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O [ Abstract 

^,„^ ■ We point out that a newly introduced recursive algorithm for lattice poly- 

^vq ■ mers has a much wider range of applicability. In particular, we apply it to the 

simulation of off-lattice polymers with Lennard-Jones potentials between non- 
\^J ' bonded monomers and either delta or harmonic potentials between bonded 

l/-^ • monomers. Our algorithm allows particularly easy calculations of the free 

g : energy, and .een« in general more efficient than oUrer existing algorithms. 
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1 Introduction 

In a number of recent papers, new Monte Carlo schemes for simulating off-lattice 
polymers have been proposed ll|]^[il- A particular aspect there was the calculation 
of chemical potentials. This is not easy in most schemes, which has lead to some 
controversy P, 0]- In this note we want to point out that a recently introduced 
recursive implementation of the enrichment method for lattice polymers [p 



can 



be easily adapted to this problem. It is both easy to implement and efficient, and 
the computation of the chemical potential is straightforward. Indeed, for polymers 
without long-range monomer- monomer interactions and without interactions with 
any solvent our method seems to be faster than all methods mentioned above. It is 
not efficient for systems at extremely low energies (when Boltzmann factors due to 
attractive potentials between individual monomers become > 10), and for systems 
with long-range forces ll^, [l3l . 

Though we shall apply our algorithm only to polymer systems, we should point 
that it is much more general. It can be applied to any equilibrium system which 
can be broken up into discrete units, labeled by an index i = 1, . . . N, and whose 
internal energy can be written as 

N 

U = Y.Ur. (1) 

1=1 

where Ui depends only on units with label i' < i. In the case of a polymer with 
potential Uij between non-bonded monomers and potential Vi between monomers i 
and i — 1, we choose of course Ui = Vi + J2j<i Uij for i > 1, while Ui = Vi is evaluated 
with xq = 0. Thus the start of polymer chain is anchored to a: = 0. 

The algorithm basically tries to build the system by assembling it unit by unit. 
To obtain the correct Boltzmann weights, the entire assembled configuration has 
to be discarded occasionally (with probabilities depending on these weights) if the 
following units do not fit. A typical example is a self-avoiding walk on a lattice, where 
the configuration has to be discarded if the next step would lead to a self intersection. 
In order to compensate for this "attrition", we use an "enrichment" method, the 
basic idea of which was already introduced more than 30 years ago [|14|. Instead 
of trying just one continuation from a partially assembled system, the intermediate 
sample is "enriched" by replicas which serve as starting points for independent 
continuations. Thus from each partly assembled system more than one continuation 
is attempted. These attempts are made irrespective of whether one of them is 
successful or not, which leads immediately to an unbiased Gibbs ensemble and which 
distinguished the method from most other proposals to overcome attrition. 

But in contrast to earlier implementations of the enrichment idea which were 
of "breadth first" type |]T5|, |T6|, |1^, [1^, our implementation is "depth first" p!9| . 



This implies completely different data structures. In particular, it means that a 
partly assembled system can reside in the fast memory of the computer even in 
very large simulations. This avoids the very time consuming data transfer needed 
in breadth-first implementations, unless the systems are very small. As discussed in 



17], simulating a system with N units one has to simulate >> N ensemble members 



simultaneously in a breadth first algorithm. This requires a storage space > (9(A^^), 
which limits severely the system sizes. But even if the entire ensemble fits into 
memory, a breadth first algorithm is slower by a factor 0{N) since adding a new 
unit takes a time 0{N) [|r^ (there is a finite probability that a new replica has to 
be created) while it needs only a time 0{1) for a depth first algorithm. Also, the 
simplest and most intuitive implementation of a depth first algorithm is via recursive 
function calls. In this case the compiler makes all the book-keeping which is fast 
but quite non-trivial in this approach. Thus our method avoids all problems which 
have made the enrichment method unpopular in the past. 

The only disadvantage of a depth first implementation is that we have to guess 
the attrition in advance. In a breadth-first approach, level i is treated only after 
all previous levels have been finished, and thus the attrition on the previous levels 
is known. If we know the amount of attrition sufficiently well from other sources, 
we do not have any problem in a depth first approach either. Otherwise, the best 
strategy is to start with small systems and a conservative estimate of the attrition, 
and to increase the system size in separate runs as the attrition gets better and 



better known 1 20 



2 Algorithm 

Our aim is to compute the partition function 

Zjv = I dxi . . . rfx^ e-^u{.,....^) ^ p ^ (kBT)-^ . (2) 

In addition, we consider also "partial partition functions" which describe the last 
N — i units in the static background field created by the first i + 1 units, 

Z^_,|,(xi ... X,) = I rfx,+i ...d^N e-^^r^^+i ''^^'■■■''^^ . (3) 

They can be written recursively, 

Z^_,+i|i„i(xi, . . . ,Xi_i) = y dxi e-^^'(^i-''') Zjv-i|i(xi, . . . ,Xi) (4) 

with Zo|7v(xi, . . . ,XAr) = 1 and Zi\f\o = Z^. The basic strategy will be to compute 
Monte Carlo estimates for the partial partition functions using this recursion relation 
(this will be done implicitly, and the explicit code needed to do it will be very 
compact). The total partition function is generated by "assembling" the last units 
first (which means just summing over suitable statistical samples), and working one's 
way back. The task is completed when finally the sample points for the first unit 
are summed over. With this strategy in mind we will in the following concentrate 
on one typical step in the recursion relation where Zjv-j|j(xi . . . Xj) is assumed to be 
known, and Zjv-j+i|i-i(xi . . .Xj_i) is to be computed. 

We assume that the potential Ui can be split into two parts, 

f/,(xi . . .X,;) = [/(°)(xi . . .X,) + Af/,(Xi . . .X,) , (5) 



with the following properties: 

(i) The partial partition functions Zj^_^,^ associated with f/^^^ can be computed 
for each background configuration (xi . . . x,) either analytically or by some other 
method which is independent of the present Monte Carlo method; 

(ii) A fast (pseudo-) random number generator exists which produces points $,k 
distributed with the density 

(notice that this is normalized due to eq.(^); 

(iii) AUi is "kind" in the sense that the integral J d^pl 'e~^^^^ converges for 
all background configurations and is not too big. We should stress that the last 
condition affects only the efficiency of the method, but is not related to any bias. In 
particular, we make no series expansion or truncation in higher powers of Af/j, or 
anything of that sort. In general, it will be sufficient if Ul has correct asymptotic 
behavior for the integral to converge, and has roughly the same shape as f/j. 

Our U^^' is similar to the "guiding field" in [0, |I^. 

Using this decomposition of Ui one shows easily that 

Z^_i+i|i_i(xi...Xi_i) = zJJ'^^i|._i(xi...Xi_i)y ci^pf^(^|xi...Xi_i) (7) 

7(0) , 

'N-i\i\ 



-^Af-iu(^l • • -Xj-ljO 



The integral over ^ can now be approximated by a sum over random points C,k 
obtained by means of the above random number generator, and we obtain 

1 M 

ZN-^+l\^-l{^l . . . Xi_i) = ^5J^+i|,_i ^lini^ — Yl (9) 

Assume that we have already an estimator for Z^-iii- Then an estimator for 
Zjv-j+i|i-i is obtained by either associating a weight 

W(e,-X, X- i) - — ^^^-^+11'-!^''^ • ■■'^'-^^ -/3A^.(xi...x,.g,) ^^^ 

with each ^k, or — and this is the method used in our approach — by replacing 
each ^k in the average by ptWi replicas of itself, each counted with weight 1 and 
labeled by an index a. Here pi is a parameter (independent of Xj and ^k) which is in 
principle arbitrary (more about it will be said below) and which controls the size of 



the sample by counterbalancing the attrition^ during the step i ^ i — 1. This gives 
then our MC estimate 

1 ^^ 

2'Ar_i+i|i_i(xi...Xi_i) ^ — ^ J2 ^^Li|i(xi...Xi_i,^fc) • (12) 

■r J k = l replicas a 

The superscript a on the partial partition function on the rhs. is to indicate that 
we use of course different sample points (xj+i . . . x^v) for each replica, even though 
the backgrounds are the same. 

The parameter pi has to be chosen carefully: large Pi will lead to samples whose 
sizes increase quickly with i, leading to excessive CPU times for large A^, while small 
Pi lead to too small samples for large i. Optimally, Pi should be chosen such that 
the sample size is roughly independent of i. Notice that this means in particular 
that M will not be large. Indeed, most numerical results quoted below are obtained 
with M = 1. Large statistics is not obtained by trying many different continuations 
of each partially built chain, but by making many independent runs. It is only for 
extremely low temperatures (not studied here) that M >> 1 should be advantageous 
since it allows a more uniform covering of configuration space. 

We should point out that the above two possibilities (of using Wi either as a 
weight or as a multiplicity of replicas) are indeed just two special cases within a 
much wider range of possible choices. They are all distinguished by a different 
balance between equidistribution of the statistical sample and equal weights put on 
all sample point. It might well be that for different purposes different variants are 
optimal, but we shall in the following discuss only the choice of uniform weights, 
corresponding to perfect importance sampling. It has the advantage that all thermal 
averages are just normal averages without any additional weight factors except for 
the weights pi. In particular, the incremental chemical potential is simply 

^, = -ksT log -^^-ksT log ^^^^ (13) 

A-i mi^ipi 

where rrii is the total number of sample point replicas at level i (i.e., the total number 
of subroutine calls at depth i in the algorithm described below). 

Technically, our method is implemented by means of a recursively called subrou- 
tine which has the level i as argument. Basically, it just chooses a random point ^, 
inserts a monomer at this point and computes its weight u'j(^), and makes in the 
average PiWi{^) calls to itself if i < A^, with new argument i + 1. After returning from 
all these subroutine calls, the monomer at ^ is removed and the subroutine is left. 
This is of course complemented by updating the statistics for whatever observable 
is to be measured. 

Finally, we have to specify what we mean by "make ... calls in average". In 
principle, we can choose any distribution for the number of calls, provided it gives 
the right average. But efficiency will depend on this choice. One possibility would 



"'^We use the word "attrition" for conformity with the use in the literature on self- avoiding walks. 
It should be noted that in our case it does not necessarily imply a depletion of chains, but can also 
imply the opposite, depending on the sign of the potential. 



be e.g. a Poissonian (this would be in spirit with the first apphcation of this method 
in a lattice model [Q), but in general this is not the best choice. As in lattice 
models, it seems that the optimal choice depends on the specific situation, and in 
extreme cases (steep potentials and low temperatures) some experimentation will 
be needed. As a rule of thumb we propose to choose the distribution such that it 



has the smallest possible variance [|TT|. Thus, if PiWi{^) = k + rj with rj G [0, 1) and 



k integer, we first make k calls and then add one more call with probability rj. 

3 Applications 

We applied this to two versions of a model where non-bonded monomers interact 
by Lennard- Jones potentials in 3-d space, 

«u = 4K^^^ - r^j^] , rij = \rij\ = |xi - Xj-| . (14) 

In the first version 0, the force between bonded monomers is harmonic, 

v^='^ (r,,,_i - 1)2 , K = 400 (15) 

for rj^j_i > 0.5. For rj^j_i < 0.5, the potential is infinite. In the second version 0, 
it is simply provided by hard rods which keep a fixed distance r^ j_i = 1. Notice 
that we did not truncate the potential at large distances, in contrast to previous 
studies ||l|, |3]- As far as we can see, the approximations involved in correcting for 
this truncation should be the only possible source for eventual disagreements with 
these works. 

In the second version, a natural choice of t/j is given by an isotropic delta 
potential^ fixing r^ j_i at 1. Thus the vectors C, were chosen randomly on the surface 
of a sphere centered at Xj_i. The corresponding partial partition functions are 

In the first version the choice is less obvious since it is not so easy to produce 
random points according to e~^^' with Vi given by eq.(|T3p. We thus took 

ui'^=v,in,^,) + kBT\n{rl_,) (16) 

For this choice the radial distribution function is a Gaussian centered at r = 1, 
pW(ri,i_i) oc rrf_ie-^^% and Z^°) = [(27r)V(fi:/?)]"/2 for K/3»lf\ 

Fig.|l] shows our data for the excess chemical potential /3/i®^ for the hard rod and 
the harmonic potential, respectively. The value f3 = 1/1.2 was chosen to compare 
our data to those of 0]. From fig. 5 of that paper we see that Pfi"^"^ = —10.4 ± 0.4 for 
a hard rod chain with N = 30 monomers. This is much smaller (in absolute value) 
than our estimate/3/i^^ = —15.35 ± 0.05. Part of this discrepancy can be explained 



^By this we mean of course not strictly a delta function but a potential which is enough singular 
to give a delta function for the equilibrium density, e~^^i — (5(ri.i_i — 1). 

■^If this condition is not fulfilled we get an additional term proportional to an error function due 
to the centering of ri^i-i around 1. 
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Figure 1: Results for the excess chemical potential at T = 1//3 = 1.2. The upper curve 
is related to the hard rod potential, while the lower one is the result for the harmonic 
potential of eq.(^). Statistical errors are smaller than the thickness of the line. 



by the fact that the Lennard- Jones potential was truncated in at r = r^ = 2.5 and 
its contribution from r > r^ was estimated analytically. This was not done in our 
simulations, where all potentials were taken into account exactly. But by performing 
the same truncation as in 0] without correcting for it at all, we estimated that even 
this is much too small (^ 9 — 10%) to explain the discrepancy. The simulations of 
[|I|] have too large statistical errors for a similarly detailed comparison. 

Averages of the squared end-to-end distance for longer chains are shown in fig. |^. 
They clearly indicate that (3 = 1/1.2 is far below the G-temperature in agreement 
with the fact that /i'^^ is negative and decreasing with A^. In fig. ^ we show our data 
for /3/i™*^ = /3(/i^ — /i^_i), which represents the free energy needed to add one more 
monomer to the chain. At the 0-point we expect /j,^^'^ to be independent on A^. 
The chains are still too short to pin down (3q exactly, but the plots unambiguously 
show that Pq is much smaller then 0.36, the value given in [0. Together with the 
data from fig. || we would say that 0.2 < /^e < 0.23. Both, the data of fig. |^ and 
that of fig. ^ were produced using the hard rod potential for bonded monomers and 
the full LJ potential for nonbonded monomers, i.e. without performing any cutoff 
at large Vij. Each data set is based on at least 10^ 'tours'. We call a tour the set of 
all (correlated) chains produced by the same initial subroutine call from the main 



routine. 
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Figure 2: Simulation results for the average squared end-to-end distance for /3 = 
0.175, . . . , 0.375 with A/3 = 0.025 and additionally for /3 = 0.4, 0.6, 1/1.2. The data clearly 
show that the chains are already collapsed at /? > 0.25. 



4 Generalizations and Outlook 



In the above applications, we did not truncate the potential at large r. Thus, insert- 
ing a new monomer takes a time 0{N). For larger N this is no longer tolerable. If 
the potential is truncated in such a case, one should also use efficient data structures 
for searching relevant neighbors [BlI]. With this we can achieve 0{1) behavior. 



Our algorithm as described above can become inefficient for two main reasons, 
but in both cases this can be cured by minor modifications: ffist of all, if the 
temperature is very low, the Boltzmann factors e~^^^' for single monomers can 
become very large. The efficiency of the method results from the fact that large 
Boltzmann factors for the entire chain are split into smaller factors for individual 
monomers. If the latter factors themselves become large and are spread over a large 
interval, then most random positions will either be immediately discarded since they 
have very small Wj, or they will have very large Wi and lead thus to very many replicas 
and correspondingly to huge fluctuations. A case where this made our method quite 



unsuccessful is the random heteropolymer model of |22 . 
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Figure 3: Simulation results for IS^j}'"'' for f3 = 0.175, . . . ,0.375 with A/3 = 0.025. One 
can see that //™'^ is independent of N for f3 ~ 0.2 to 0.225. So we can expect the 0-point 
somewhere between these two values. 



The simplest way out of this dilemma consists in choosing M >> 1. In this 
way we will have even more points which are essentially useless since they have very 
small Wi. But the regions in state space with high probability will be sampled more 
evenly, and this should be more important. Another popular remedy consists in 
making the bulk of the simulations at higher temperatures, and quenching down to 



the desired temperatures in regular time intervals [p!q] . 

The other case where the above algorithm gives poor results is provided by 
systems in which originally favored configurations lead finally to dead ends, while 
originally unfavored configurations become more important as the growth of the 
system continues. Two specific examples are dense polymer systems in 2-d M and 



polymers with long range repulsive forces |T^, [T^ . 

Assume we want to grow a long self avoiding walk in a finite 2-d region such that 
it fills a large fraction of the area. As we place the first monomers, all configurations 
are equally likely. But those which leave large closed voids all are dead ends since 
the walk later on cannot penetrate into the voids. Obviously it would be much 
better in this case to bias the walks against large voids from the very beginning. 

The situation is similar for a polymer with repulsive Coulomb potentials. There, 
unless the force is very strong, the effect of the repulsion is not too strong for end 
monomers, and hence new monomers are added without a strong radial bias. But the 
stretching force on a monomer deep inside the chain is much bigger (since all forces 
essentially add up), and the configuration is much more stretched inside the chain. 
Thus, when an existing chain is to be prolongated, most of the existing configurations 



have to be discarded, in order to be replaced by stretched configurations which at 
first (when they are assembled) are very improbable. 

As we said, a way out of this problem is to use biased walks. This means that 
the number of replicas is not strictly proportional to Wi but is larger in those regions 
which we suspect to become more important later. Of course it means also that we 
have to replace eq.(|^) by a weighted sum. It is e.g. known that the Rosenbluth 
trick 1^ leads to a bias towards more compact SAW's. We found indeed that our 
method with Rosenbluth weights instead of uniform weights was more efficient in 
giving SAW's which fill a square with periodic boundary conditions, but it leads to 
much larger statistical fiuctuations at low density. We should point out that this 
possibility of biasing is independent of the choice of U^^' , something which seems to 
have been missed in [|16], 0, p!8[| . 

Finally, we tried our method also for non-polymeric systems. For instance, we 
simulated the 2-d Ising model with spins numbered in the same way as they would be 
e.g. stored in a FORTRAN array. Though the method worked decently, it could not 
compete either with cluster fiip methods (due to their much more efficient moves) 
nor with conventional Monte Carlo schemes which can be made very efficient by 
vectorization and multi-spin coding. We should mention that the possible applica- 
tion spin models, and to the Ising model in particular, was also pointed out in [Q 
in the context of the breadth-first approach. 
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Unfortunately, this happened only after we had used our method in several papers 
(and after the present paper was essentially finished). We apologize to the above 
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